//==============================================================================
// BLASTWrapper.cc
//------------------------------------------------------------------------------
// 
//------------------------------------------------------------------------------
// $Id: $
//------------------------------------------------------------------------------
// $Log: $
//   
//==============================================================================
//.......1.........2.........3.........4.........5.........6.........7.........8
//

#include "BLASTWrapper.hh"
#include "LogStream.hh"

//extern LogStream ls;

//string gHTMLFeatureNames[]={"Score","a name =","Query=","END"};
//string gFeatureNames[]={"Score",">","Query=","END"};

string gResultRoot = "/tmp/kjd.BLASTWrapper.out";
string gQueryRoot = "/tmp/kjd.BLASTWrapper.query";


//---------------------------------------------------------------------------
// * BLASTWrapper
//--------------------------------------------------------------------------- 
//
BLASTWrapper::BLASTWrapper()
{
  mType = "blastn ";
  mView = " -m0 ";
  mOutFileName = "";
  mQueryFileName = "";
  mParseType = "m0";

  InitFeatureNames();
}

//---------------------------------------------------------------------------
// * CreateTemporaryFileName
//--------------------------------------------------------------------------- 
// 
//
string
BLASTWrapper::CreateTemporaryFileName(char* RootName)
{
  string temp = RootName;
  return(CreateTemporaryFileName(temp));
}

//---------------------------------------------------------------------------
// * CreateTemporaryFileName
//--------------------------------------------------------------------------- 
// Adds a random suffix to a file name.  The goal is to minimize the 
// chance of conflicts with existing files from this code running 
// simultaneously as a different process on the same machine.  This 
// is still no guarantee, but it's a simple scheme that will keep 
// the number of conflicts low. 
// 
string
BLASTWrapper::CreateTemporaryFileName(string &RootName)
{
  // Give temporary file a random suffix to minimize chance of 
  // collisions.  KJD can we do better???
  RandomNumberGenerator RNG;
  RNG.SeedWithClock();
  long rndSuffix = (long) RNG.Uniform(0,RAND_MAX);
  char buf[80];
  sprintf(buf,"%ld",rndSuffix);
  
  string OutName = RootName+buf;
  return(OutName);
}

//---------------------------------------------------------------------------
// * RunBLASTOnFile
//--------------------------------------------------------------------------- 
// Constructs a BLAST command and executes it. 
// 
void 
BLASTWrapper::RunBLASTOnFile()
{
  long StartTime;
  StartTiming(StartTime);
  
  string BLASTCommand = 
	 "blastall -p "+mType+" -d "+mBLASTDBName+
	 " -i "+mQueryFileName+mView+" > " +mOutFileName;
  
   system(BLASTCommand.c_str());
	mElapsedTime = EndTiming(StartTime);
}

//---------------------------------------------------------------------------
// * ParseBLASTOutput
//--------------------------------------------------------------------------- 
// Coordinates parsing BLAST output file into a vector of results. 
//
void
BLASTWrapper::ParseBLASTOutput(BLASTSearchResultVec &ResultsVec)
{
  if (mParseType == "m0"){
	 ParseView0(mOutFileName,ResultsVec);
  }else{
	 if (mParseType == "m0html"){
		ParseView0HTML(mOutFileName,ResultsVec);
	 }else{
		cout<<"Unknown parse type specified:"<<mParseType<<endl;
	 }
  }
}

//---------------------------------------------------------------------------
// * MapSamplesToConcatenated
//--------------------------------------------------------------------------- 
// Takes the result vector from a BLAST on concatenated samples (done for 
// BLAST speed) and properly maps them back onto the original samples.  
// Original sample names are saved in the result vec.  
//
void
BLASTWrapper::MapSamplesToConcatenated(BLASTSearchResultVec &ResultsVec)
{
  for(int r = 0;r < ResultsVec.size();r++){
	 long QLoc = ResultsVec[r].mQueryLocation;
	 int SampleNumber = QLoc/mSampleSize;
	 
	 long SEnd0 = ((SampleNumber+1)*mSampleSize)-1;
	 long SEnd1 = ((SampleNumber+2)*mSampleSize)-1;
	 long QEnd = ResultsVec[r].mQueryLocation+ResultsVec[r].mIdentitiesDen-1;

	 //cout<<"r="<<r<<" QLoc="<<QLoc<<" SampleNumber="<<SampleNumber<<endl;

	 // Sanity check
	 if (SampleNumber >= mpSamplesVec->size()){
		cout<<"ERROR SampleNumber="<<SampleNumber<<" VecSize=";
		cout<<mpSamplesVec->size()<<endl;
	 }else{

		// Handle the case where a sample hit straddles a sample boundary 
		// in the concatenated sequence. 
		if (SEnd0 < QEnd){
		  // If it doesn't straddle more than two samples, then create two 
		  // sequences out of it.  Since it is almost always the case that
		  // one sequence will be large and the other very small, we only 
		  // save the larger one. 
		  if (QEnd <= SEnd1){
			 long SLen1 = SEnd0 - QLoc + 1;
			 long SLen2 = QEnd - SEnd0;

			 //ResultsVec[r].mbSplitConstructed = true;

			 if (SLen1 > SLen2){
				ResultsVec[r].mQueryName = (*mpSamplesVec)[SampleNumber].mDescription;
				ResultsVec[r].mIdentitiesNum = SLen1; // KJD !!! Not necessarialy true
				ResultsVec[r].mIdentitiesDen = SLen1; 
			 }else{
				ResultsVec[r].mQueryName = 
				  (*mpSamplesVec)[SampleNumber+1].mDescription;
				ResultsVec[r].mQueryLocation = SEnd0+1;
				ResultsVec[r].mIdentitiesNum = SLen2; // KJD !!! Not necessarialy true
				ResultsVec[r].mIdentitiesDen = SLen2; 
			 }
			 //cout<<"Split: ";
			 //mResultsVec[r].Print(cout);
		  }else{
			 //cout<<"Double split:"<<endl;
		  }
		}else{
		  // Hit only encompases one sample, so simply need to record the 
		  // name of the original source sample.
		  //cout<<"Contained hit:";
		  ResultsVec[r].mQueryName = (*mpSamplesVec)[SampleNumber].mDescription;
		  //ResultsVec[r].mbSplitConstructed = false;
		  //mResultsVec[r].Print(cout);
		}
	 }
  }
}



//  for(int r = 0;r < mResultsVec.size();r++){
//	 int SampleNumber = mResultsVec[r].mQueryLocation/mSampleSize;
//	 mResultsVec[r].mQueryName = (*mpSamplesVec)[SampleNumber].mDescription;
//  }


 
//---------------------------------------------------
// * ConcatenateSamples
//---------------------------------------------------
// Takes a vector of sequences and constructs a single sequence
// which is the concatenation of those sequences.  BLAST incurs a 
// penalty for each query sequence, so a substantial performance 
// improvement is obtained by concatenating a large sample set before 
// calling BLAST.  
//
void 
BLASTWrapper::ConcatenateSamples(Sequence &target,SequenceVector &Samples)
{
  Sequence n;
  
  // create a block of n's to pad between samples 
  for(int j = 0;j<20;j++){
   n.push_back('n');
  }

  for(int i = 0; i< Samples.size();i++){
	 target+=Samples[i];
	 //target+=n;  // KJD Bug test, insert block of n's
  }
  target.mDescription = "Concatenated query.";
}


//---------------------------------------------------
// * StripReadNameSuffixes
//---------------------------------------------------
//
void 
BLASTWrapper::StripReadNameSuffixes()
{
  for(int r = 0;r < mResultsVec.size();r++){
	 int idx = mResultsVec[r].mQueryName.find(".");
	 if (idx >=0){
		string temp = mResultsVec[r].mQueryName.substr(0,idx);
		cout<<mResultsVec[r].mQueryName<<"->"<<temp<<endl;
		mResultsVec[r].mQueryName = temp;
	 }
  }
}

//---------------------------------------------------
// * StripAccessionNumber
//---------------------------------------------------
//
void 
BLASTWrapper::StripAccessionNumber()
{
  for(int r = 0;r < mResultsVec.size();r++){
	 if (mResultsVec[r].mQueryName.size()>=8){
		string temp = mResultsVec[r].mQueryName.substr(1,8);
		cout<<mResultsVec[r].mQueryName<<"->"<<temp<<endl;
		mResultsVec[r].mQueryName = temp;
	 }
  }
}


//---------------------------------------------------
// * SearchConcatenatedSamples
//---------------------------------------------------
// Coordinates the steps in searching a BLAST database with a given
// set of samples, concatenated for speed. 
//
//
void 
BLASTWrapper::SearchConcatenatedSamples(SequenceVector &samplesVec,
					string DBName,BLASTSearchResultVec &ResultsVec)
{
  mpSamplesVec = &samplesVec;
  mBLASTDBName = DBName;
  mQueryFileName = CreateTemporaryFileName(gQueryRoot);
  mOutFileName = CreateTemporaryFileName(gResultRoot);

  // Concatenate samples into a single sequence
  // KJD What a horrible kludge to write a single sequence
  SequenceVector ConcatenatedSamples;
  ConcatenatedSamples.GrowVec();
  ConcatenateSamples(ConcatenatedSamples[0],samplesVec);

  // Save the size of the samples for later parsing. 
  mSampleSize = samplesVec[0].size(); 

  FASTAFile FAFile;
  FAFile.Write(mQueryFileName,ConcatenatedSamples);

  RunBLASTOnFile();

  // Remove the temporary sample file...
  string rmCommand = "rm -f "+mQueryFileName;
  system(rmCommand.c_str());

  ParseBLASTOutput(ResultsVec);

  // Map the concatenated sequence back to the samples. 
  MapSamplesToConcatenated(ResultsVec);

  if (mbViewByRead){
	 StripReadNameSuffixes();
  }

  // KJD Kludge
  StripAccessionNumber();


  // Remove temporary BLAST Output file created in RunBLASTOnFile...
  rmCommand = "rm -f "+mOutFileName;
  system(rmCommand.c_str());

 // Extract subject and query name information
  //cout<<"get subject list: "<<flush;
  ResultsVec.GetSubjectList(ResultsVec.mSubjects);
  //  cout<<ResultsVec.mSubjects.size()<<"..."<<flush;
  //cout<<"get query list: "<<flush;
  ResultsVec.GetQueryList(ResultsVec.mQueries);
  //cout<<ResultsVec.mQueries.size()<<"..."<<flush;
  //cout<<"sort..."<<flush;
  sort(ResultsVec.mSubjects.begin(),ResultsVec.mSubjects.end());
  sort(ResultsVec.mQueries.begin(),ResultsVec.mQueries.end());

}

//---------------------------------------------------------------------------
// * Search
//---------------------------------------------------------------------------
// BLAST against a vector of sequences. 
//
void
BLASTWrapper::Search(SequenceVector &seqVec,string DBName,
		     BLASTSearchResultVec &ResultsVec)
{
  mQueryFileName = CreateTemporaryFileName(gQueryRoot);
  mOutFileName = CreateTemporaryFileName(gResultRoot);
  mBLASTDBName = DBName;

  cout<<"Creating file "<<mQueryFileName.c_str()<<" with "<<seqVec.size();
  cout<<" sequences for BLAST."<<endl;

  FASTAFile FAFile;
  FAFile.Write(mQueryFileName,seqVec);


  cout<<"BLAST...."<<flush;
  RunBLASTOnFile();
  cout<<"run complete. "<<endl;

   // Remove the temporary sample file...
  string rmCommand = "rm -f "+mQueryFileName;
  system(rmCommand.c_str());

  cout<<"Parsing BLAST output..."<<flush;
  ParseBLASTOutput(ResultsVec);
  cout<<"done... "<<ResultsVec.size()<<" results found."<<endl;

  // KJD Stripping Read Name suffixes should probably happen elsewhere
  //if (mbViewByRead){
  //	 ls<<"Strip read name suffixes (to view by read)..."<<flush;
  //	 StripReadNameSuffixes();
  //	 ls<<"done."<<endl;
  //}

  // Remove temporary BLAST Output file created in RunBLASTOnFile...
  rmCommand = "rm -f "+mOutFileName;
  system(rmCommand.c_str());

  // Extract subject and query name information
  //cout<<"get subject list: "<<flush;
  ResultsVec.GetSubjectList(ResultsVec.mSubjects);
  //cout<<ResultsVec.mSubjects.size()<<"..."<<flush;
  //cout<<"get query list: "<<flush;
  ResultsVec.GetQueryList(ResultsVec.mQueries);
  //cout<<ResultsVec.mQueries.size()<<"..."<<flush;
  //cout<<"sort..."<<flush;
  sort(ResultsVec.mSubjects.begin(),ResultsVec.mSubjects.end());
  sort(ResultsVec.mQueries.begin(),ResultsVec.mQueries.end());

}

//---------------------------------------------------------------------------
// * Search
//---------------------------------------------------------------------------
// Perform BLAST search against already existing file. 
// 
void
BLASTWrapper::Search(string TargetName,string DBName,
							BLASTSearchResultVec &ResultsVec)
{
  mQueryFileName = TargetName;
  mOutFileName = CreateTemporaryFileName(gResultRoot);
  mBLASTDBName = DBName;

  RunBLASTOnFile();

  ParseBLASTOutput(ResultsVec);

  if (mbViewByRead){
	 StripReadNameSuffixes();
  }

  // Remove temporary BLAST Output file created in RunBLASTOnFile...
  string rmCommand = "rm -f "+mOutFileName;
  system(rmCommand.c_str());

  // Extract subject and query name information
  ResultsVec.GetSubjectList(ResultsVec.mSubjects);
  ResultsVec.GetQueryList(ResultsVec.mQueries);
  sort(ResultsVec.mSubjects.begin(),ResultsVec.mSubjects.end());
  sort(ResultsVec.mQueries.begin(),ResultsVec.mQueries.end());
}

//---------------------------------------------------------------------------
// * InitKeywordList
//---------------------------------------------------------------------------
// 
void
BLASTWrapper::InitKeywordList2(string Names[])
{
  cout<<"New version of InitKeywordList"<<endl;

  int idx=0;
  string temp;
  if (mKeywordList.size()>0){
  	 mKeywordList.erase(mKeywordList.begin(),mKeywordList.end());
  }

  cout<<"Starting first list"<<endl;
  while(Names[idx] != "END"){
	 cout<<Names[idx]<<endl;
	 idx++;
  }
  cout<<"Done with first loop."<<endl;
  idx = 0;
  while(Names[idx] != "END"){
	 temp = Names[idx];  // for debugger observability
	 cout<<temp<<endl;
	 mKeywordList.push_back(temp);
	 idx++;
  }
  cout<<"Done with loop"<<endl;
}

//---------------------------------------------------------------------------
// * InitFeatureNames
//---------------------------------------------------------------------------
// 
void 
BLASTWrapper::InitFeatureNames()
{
  if (mKeywordList.size()>0){
  	 mKeywordList.erase(mKeywordList.begin(),mKeywordList.end());
  }  

  mKeywordList.push_back("Score");
  mKeywordList.push_back(">");
  mKeywordList.push_back("Query=");
}

//---------------------------------------------------------------------------
// * InitHTMLFeatureNames
//---------------------------------------------------------------------------
// 
void 
BLASTWrapper::InitHTMLFeatureNames()
{
  if (mKeywordList.size()>0){
  	 mKeywordList.erase(mKeywordList.begin(),mKeywordList.end());
  }  

  mKeywordList.push_back("Score");
  mKeywordList.push_back("a name =");
  mKeywordList.push_back("Query=");
}


//---------------------------------------------------------------------------
// * ParseView0
//---------------------------------------------------------------------------
// Parses the -m0 view of BLAST output into a vector of hit information.  
//
void
BLASTWrapper::ParseView0(string FileName,BLASTSearchResultVec &ResultsVec)
{
  
  BLASTSearchResult tempResult;
  string tempLine;
  int offset=0;
  int idx;
  
  InitFeatureNames();

  ifstream inFile;
  inFile.open(FileName.c_str(),ios::in);
  if (inFile.fail()){
	 cout<<FileName<<" file open failed (Read Multiple FASTA)"<<endl;
	 return;
  }

  GetNextLineContaining("Query=",inFile,tempLine);
  if (tempLine =="") return;   // empty file.
  GetRemainderAfter("Query=",tempLine,tempResult.mQueryName);

  GetNextLineContaining("Database:",inFile,tempLine);
  GetRemainderAfter("Database:",tempLine,ResultsVec.mDatabaseName);
  
  // Occurence of a > marks the name of a subject file...
  GetNextLineContaining(">",inFile,tempLine);
  GetRemainderAfter(">",tempLine,tempResult.mSubjectName);

  GetNextLineContaining("Length = ",inFile,tempLine);
  offset = 0;
  tempResult.mSubjectLength = (long) GetNextDouble(tempLine,offset,offset);

  do{
	 GetNextLineContaining("Score",inFile,tempLine);
	 int KeyWordIdx=0;
	 do{
		offset=0;
		double temp = GetNextDouble(tempLine,0,offset); // Raw score
		tempResult.mScoreBits = GetNextDouble(tempLine,offset,offset); 
		tempResult.mExpect=GetNextDouble(tempLine,offset,offset);      

		GetNextLineContaining("Identities",inFile,tempLine);
		
		tempResult.mIdentitiesNum = GetNextInteger(tempLine,0,offset);
		tempResult.mIdentitiesDen = GetNextInteger(tempLine,offset,offset);
		
		GetNextLineContaining("Query:",inFile,tempLine);
		tempResult.mQueryLocation = GetNextInteger(tempLine,0,offset); 
		
		GetNextLineContaining("Sbjct:",inFile,tempLine);
		tempResult.mSubjectLocation = GetNextInteger(tempLine,0,offset);

		// Add results to the list
		// KJD: This crude parser grabs some random lines with "Score" on
		// the line.  Those lines return mScoreBits==0, which is not a 
		// valid BLAST score, so I simply omit those here. 
		if (tempResult.mScoreBits != 0){
		  ResultsVec.push_back(tempResult);
		}
	 
		// Find the next keyword, which might be ">", "Score", or "Query="
		KeyWordIdx=GetNextLineContaining(mKeywordList,inFile,tempLine);
		
		// while == Score, i.e. just more hits for same subject
	 }while((mKeywordList[KeyWordIdx]=="Score") && (!inFile.eof()));
	 
	 // Found a new subject
	 if (mKeywordList[KeyWordIdx] == ">"){
		GetRemainderAfter(">",tempLine,tempResult.mSubjectName);
		GetNextLineContaining("Length = ",inFile,tempLine);
		offset = 0;
		tempResult.mSubjectLength = (long) GetNextDouble(tempLine,offset,offset);
	 }
   
	 if (mKeywordList[KeyWordIdx]=="Query="){		
		GetRemainderAfter("Query=",tempLine,tempResult.mQueryName);
		
		GetNextLineContaining("Database:",inFile,tempLine);
		GetRemainderAfter("Database:",tempLine,ResultsVec.mDatabaseName);
		
		GetNextLineContaining(">",inFile,tempLine);
		GetRemainderAfter(">",tempLine,tempResult.mSubjectName);

		GetNextLineContaining("Length = ",inFile,tempLine);
		offset = 0;
		tempResult.mSubjectLength = (long) GetNextDouble(tempLine,offset,offset);
	 }	  
  }while(!inFile.eof());

  inFile.close();
}


//---------------------------------------------------------------------------
// * ParseView0HTML
//---------------------------------------------------------------------------
// Parses the -m0 view of BLAST output returned by Harley's ncbi blast 
// search script into a vector of hit information.   These files are 
// htmlized and not based on fa files, so many of the feature identifiers
// are different than for direct blast output. 
// 
// KJD: This version of the parser does not parse out subject length.
//
void
BLASTWrapper::ParseView0HTML(string FileName,BLASTSearchResultVec &ResultsVec)
{
  
  BLASTSearchResult tempResult;
  string tempLine;
  int idx;

  InitHTMLFeatureNames();
  
  ifstream inFile;
  inFile.open(FileName.c_str(),ios::in);
  if (inFile.fail()){
	 cout<<FileName<<" file open failed (Read Multiple FASTA)"<<endl;
	 return;
  }

  GetNextLineContaining("Query=",inFile,tempLine);
  if (tempLine =="") return;   // empty file.

  GetRemainderAfter("Query=</b>",tempLine,tempResult.mQueryName);

  GetNextLineContaining("Database:",inFile,tempLine);
  GetRemainderAfter("Database:</b>",tempLine,ResultsVec.mDatabaseName);
  
  // Occurence of "a name =" marks the name of a subject file, which 
  // occurs two lines later. 
  GetNextLineContaining("a name =",inFile,tempLine);
  getline(inFile,tempLine);  // skip line
  getline(inFile,tempResult.mSubjectName); // line containing the subject name
  
  do{
	 GetNextLineContaining("Score",inFile,tempLine);	   
	 int KeyWordIdx=0;
	 do{

		int offset=0;
		double temp = GetNextDouble(tempLine,0,offset); // Raw score
		tempResult.mScoreBits = GetNextDouble(tempLine,offset,offset);
		tempResult.mExpect=GetNextDouble(tempLine,offset,offset);
	 
		GetNextLineContaining("Identities",inFile,tempLine);
		
		tempResult.mIdentitiesNum = GetNextInteger(tempLine,0,offset);
		tempResult.mIdentitiesDen = GetNextInteger(tempLine,offset,offset);
		
		GetNextLineContaining("Query:",inFile,tempLine);	 
		tempResult.mQueryLocation = GetNextInteger(tempLine,0,offset);  
		
		GetNextLineContaining("Sbjct:",inFile,tempLine);
		tempResult.mSubjectLocation = GetNextInteger(tempLine,0,offset);

		// Add results to the list
		// KJD: This crude parser grabs some random lines with "Score" on
		// the line.  Those lines return mScoreBits==0, which is not a 
		// valid BLAST score, so I simply omit those here. 
		if (tempResult.mScoreBits != 0){
		  ResultsVec.push_back(tempResult);
		}
	 
		// Find the next keyword, which might be "a name=", "Score", or "Query="
		KeyWordIdx=GetNextLineContaining(mKeywordList,inFile,tempLine);
		
		// while == Score, i.e. just more hits for same subject
	 }while((mKeywordList[KeyWordIdx]=="Score") && (!inFile.eof()));
	 
	 // Found a new subject
	 if (mKeywordList[KeyWordIdx] == "a name ="){
		getline(inFile,tempLine);  // skip line
		// line containing the subject name
		getline(inFile,tempResult.mSubjectName); 
	 }
   
	 if (mKeywordList[KeyWordIdx]=="Query="){		
		GetRemainderAfter("Query=",tempLine,tempResult.mQueryName);
		
		GetNextLineContaining("Database:",inFile,tempLine);
		GetRemainderAfter("Database:",tempLine,ResultsVec.mDatabaseName);
		
		GetNextLineContaining("a name =",inFile,tempLine);
		getline(inFile,tempLine);  // skip line
		// line containing the subject name
		getline(inFile,tempResult.mSubjectName); 
	 }	  
  }while(!inFile.eof());

  inFile.close();
}












